Algebraic perturbation theory for dense liquids with discrete potentials 
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A simple theory for the leading-order correction gi(r) to the structure of a hard-sphere liquid with 
discrete (e.g. square-well) potential perturbations is proposed. The theory makes use of a general 
approximation that effectively eliminates four-particle correlations from gi(r) with good accuracy 
at high densities. For the particular case of discrete perturbations, the remaining three-particle 
correlations can be modeled with a simple volume-exclusion argument, resulting in an algebraic and 
surprisingly accurate expression for gi(r). The structure of a discrete "core-softened" model for 
liquids with anomalous thermodynamic properties is reproduced as an application. 
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Although in principle the properties of complex phys- 
ical systems can be computed to arbitrary accuracy 
through computer simulations, perturbation theory re- 
mains instrumental across quantum, classical and statis- 
tical mechanics. Its utility is apparent whenever one is 
interested in getting quantitative insights into the be- 
havior of a system under different conditions, a task that 
would otherwise require several - possibly costly - simu- 
lations spanning the desired parameters. The underlying 
assumption is that there exists a well-characterized refer- 
ence system, the properties of which are relatively simi- 
lar to that of the system of interest, so that discrepancies 
between the two can be quantified through a small pa- 
rameter expansion. The use of a perturbative approach 
is then justified, as is often the case, if the leading-order 
terms in this expansion can be computed with greater 
ease than those required by the exact approach. 

When applied to the structure of simple atomic liq- 
uids, in which the interaction potential is comprised of 
a small, smooth and short-ranged attractive component 
outside a sharply repulsive core, liquid-state perturbation 
theory [l[ takes on a especially simple form: According 
to the Weeks- Chandler- Andersen theory (see e.g. 0]), 
the structure of such systems is largely dominated by the 
repulsive part of the potential, and particularly at high 
densities it can be approximated by the structure go(r) 
of the purely repulsive fluid alone with excellent accuracy 
Q; in this rather peculiar case, perturbative corrections 
to the reference structure are hardly necessary 0, [B| . 

For more complex liquid models with radial interac- 
tions other than such specific attractions, however, one 
cannot expect to always find simple reference systems 
with this property. Interest in such models stems in par- 
ticular from recent efforts to model jf| 0, @| and explain 
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through the existence of multiple length scales in the 
interaction potential, which introduce additional repul- 
sive and attractive interactions outside the atomic core. 
In these and similar cases, perturbation theory needs to 
account for terms beyond the reference (repulsive-core) 
structure go( r )- 



Unfortunately, the leading-order correction gi(r) to 
the radial distribution function of a liquid due to an arbi- 
trary perturbing potential is in general a function of two-, 
three_j and four-particle correlations of the reference fluid 
(cf. Consequently, over the years several approx- 

imation schemes have been proposed to simplify this cor- 
rection [3 [H, [H, [HI H, Ell • Such efforts have focused 
either on the superposition approximation, which elimi- 
nates high-order correlation functions but is known to fail 
at moderately high densities 14], or on the application of 
the integral equation method with different approximate 
closures [3 [H, [l6, F7, 3, III, which typically requires 
the numerical solution of a nonlinear integral equation 
for gi(r) at each condition of interest. In the particular 
case of square well potentials, the latter has been applied 
with varying degrees of success, depending on the density 
and the range of the well 3, 3 3] ■ 

The present paper introduces simple approximations 
that result in algebraic corrections to go(r) (cf. Eq. ([5])) 
when the perturbing potential is a discrete function of r - 
i.e. comprised of "steps" - thereby avoiding the numer- 
ical computations required by previous methods while 
achieving competitive accuracies at high densities. The 
th eory is illustrated for liquids with one [3 | and two 
0, d(| steps in the interaction potential. 

To begin the derivation, let XU(r N ) be the TV- 
body perturbation potential of a homogeneous, single- 
component fluid. The two-particle density can then be 
expanded in powers of A, and a convenient expression 

(2) 

for the leading-order correction Xp\ (ri,r 2 ) to the ref- 
erence density p Q J (ri,r 2 ) can be found by directly A- 
differentiating the definition of the two-particle density 
in the canonical ensemble, yielding 



(2 
Pi 



( ri ,r 2 ) = -rf(r x ,r 2 ) \{U)^ 2) - (f/) l . (la) 



Here /3 = 1/fcT is the usual temperature parameter, (-)o 
is the ensemble average with respect to the reference 
system, and 2 ^ is the same average with the con- 
straint that particles 1 and 2 are held fixed at ri and 
r 2 , respectively. As it stands, it can be straightforwardly 
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shown that this expression is fully equivalent to the "dis- 
of Barker and Henderson 



14] 



crete representation ' 
the limit of infinitely thin shells, or to the more tradi- 
tional result in terms of many-particle densities 0, E3] > of 
course, use of either expression in analytic form presents 
a formidable challenge, as very little is known about the 
required many-particle correlation functions. 

Let us now decompose U in terms of one-particle 
energies U(r N ) — ^J2iLi u i( rN )> where Ui(r N ) = 
Y^j 7 ki u ( T ii T j) 18 the perturbation energy of particle i 
with the remaining ones, and u(rj,rj) = u(|rj — rj|) is 
the (spherically symmetric) pairwise perturbation energy 
between particles i and j; Eq. (|la[) can thus be written 
as 



Pi 2) (r 1; r 2 ) 



-0Po 2) (ri,r 2 ) (wi)o 1:2) - (ui) 
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where use was made of the indistinguishability of the par- 
ticles wherever convenient. Consider now the last term 
in parentheses. Since the inhomogeneities caused by par- 
ticles 1 and 2 occupy volumes of molecular size only, one 

(1 2) 

expects the difference (1*3)0 ' — ( u 3)o to scale down ex- 
tensively with the total volume of the fluid. The first 
approximation to be invoked in this paper consists in as- 
suming that the resulting (in principle intensive) contri- 
bution to Eq. (fib)) due to the average of 113 is vanishingly 
small in comparison to that due to itx, so that 

^ 2) (n,r 2 ) a -^(n,ii) [(u,)^ - ( Ul ) ] . (2a) 

This approximation effectively eliminates four-body cor- 
relations from Eq. (jib)) . and thus can be expressed in 
terms of pair and triplet particle densities only: 



Pi 2) (ri,r 2 ) « -/3p 2) (ri,r 2 ) <{ u(ri,r 2 ) 



+ /dr 3 tt(ri,r 3 ) 



P 3) (ri,r 2 ,r 3 ) _ p^'{n,r 3 ) 

Po 2) (n. r 2) P 



or, in terms of distribution functions, 



g x (r) w -f3go(r) 



u(r) + p dr' it(r') Ag (r'|r) 



(2b) 



(2c) 



where p — N/V is the particle density, and Ag (r'\r) = 
5o( r 'l r ) — 5o( r ') is the change in the local structure at r' 
from a given fluid particle due to the presence of another 
particle fixed at r from it. Note that Eq. ([2]) would be 
exact had the perturbation potential acted on the total 
energy of particles 1 and 2 only, i.e. XU = X(u% + u 2 )/2, 
as for example in the theory of hydrophobicity, where 
particles 1 and 2 would play the role of solutes, and 
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FIG. 1: Radial distribution function for a hard-sphere fluid 
with a single square- well of depth e/kT = 1 extending from 
rja — 1 to r/a = 1.5, corrected according to the various 
expressions for gi(r) in the paper. The particle density is 
pa 3 = 0.85. The shaded curve is the "exact" g(r) obtained by 
direct Monte Carlo simulation. The exact lst-order correction 
is computed with Eq. (|lajl , the solvation approximation with 
Eq. (|2ap . and the excluded-volume theory with Eqs. ©-(fj]|. 



one would like to compute corrections due to solute- 
solvent and solute-solute interactions to solute-solute po- 
tentials of mean force 2l|. With this analogy in mind, 
the approximation in Eq. |(5J) is essentially equivalent 
to neglecting structural rearrangements due to "solvent- 
solvent" perturbations (i.e. interactions between parti- 
cles 3,4, . . . , N); consequently, it is expected to be more 
accurate at high densities, wherein the relative arrange- 
ment of such particles is primarily determined by hard- 
sphere packing alone. These observations motivate refer- 
ring to Eq. ^ as the solvation approximation to <?i(r). 

The validity of Eq. ([2]) has been tested for both dis- 
crete and continuous potentials typically used in liquid- 
stateperturbation theory - such as the square-well liq- 



uid |l4j, and the Barker-Henderson and Weeks-Chandler- 



Andersen partitioning schemes for Lennard-Jones poten- 
tials Q - as well as for the discrete "core-softened" model 
of Ref. As illustrated in Figs. Q] and [5] for the dis- 
crete models of present interest, the solvation approxi- 
mation provides good quantitative results for both one- 
and two-step potentials at high densities (pa 3 = 0.85); in 
agreement with the above discussion, however, the qual- 
ity of the approximation decreases at lower densities (cf. 
Fig. U and discussion on thermodynamics below). Al- 
though not central to the algebraic development below, 
it is also worth mentioning that the sampling convergence 
of Eq. (f2aj) seems to be much faster than that of the exact 
expression, Eq. (jTaJ) : it would be interesting to explore 
such issues in the context of the unbiased estimators of 
Ref. [13. 

For the particular case of discrete potentials, the sol- 
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FIG. 2: Radial distribution function for the double square- 
well potential of Ref. [j| at per 3 = 0.85. Legends as in Fig- 
ure [1] The step depths are Ei/kT = 0.5 and £2/ 'kT = 1, cor- 
responding to the steps ending at n/cr = 1.4 and r^jo = 1.7, 
respectively (inset). The excluded- volume theory is given by 
Eq. © with M = 2. 



vation approximation is considerably simplified. Con- 
sider for example the square- well perturbation to a hard- 
sphere fluid where u(r) = —e for r < r* , and u(r) = 
otherwise, with r* > a. Inserting this expression in 



Eq. © yields 



gi (r) = -(3g (r) [u(r) 



eAn*(r)} 



(3) 



where AnJ(r) is the net change in the average number of 
particles inside a shell of radius r* surrounding a given 
particle, due to the presence of another particle fixed 
at a distance r from it. The problem of computing the 
leading-order correction to go (r) for such systems thus 
reduces to that of quantifying population changes in the 
coordination shell defined by the radius r* due to an "in- 
truder" particle at r in the reference fluid. In the follow- 
ing, a simple volume- exclusion theory for Auq (r) will be 
presented. 

For the sake of clarity, let us focus on the single square- 
well case with r* = 1.5a, and let us monitor the average 
number of particles within the r*-shell of particle 1 as 
particle 2 is brought to a distance r from it (see Fig. [3]). 
To a first approximation, appreciable changes in the av- 
erage number of particles inside the r*-shell of particle 1 
should occur whenever particle 2 occupies a volume that 
would otherwise be occupied by the particles inside that 
shell. At sufficiently high densities, wherein the particles 
are tightly packed around each other, this should hap- 
pen progressively as the body of particle 2 infiltrates the 
r*-shell; a simple model that captures this idea is thus 

v*(r) 



AnS(r) 



(4) 



where v* (r) is the two-sphere intersection volume shown 
in Fig. [3J and v a = ira 3 /6 is the volume of a sphere 



FIG. 3: Two-dimensional representation of the excluded- 
volume model [Eq. (j4|] for the single square-well potential 
with r* = 1.5cr. Solid circles represent the hard spheres of 
the reference system with diameter er, while the larger dashed 
circle represents the outer radius r* of the perturbation well. 
The present model assumes that in the reference hard-sphere 
system, the average number of particles expelled from the 
r*-shell of particle 1 due to the presence of particle 2 is pro- 
portional to the volume v*(r) of particle 2 inside the shell 
(shaded area). The normalization is such that exactly one 
particle is expelled from the r*-shell when particle 2 is fully 
inside the well region delimited by the dashed circle. 



of diameter a. Note that this expression is normalized 
so that Ano(cr) = —1, i.e. exactly one particle is ex- 
cluded from the r*-shell when particle 2 is fully inside 
the r < r* region. Since a polynomial expression for the 
intersection volume v*(r) of the spheres can be obtained 
by simple quadrature, Eqs. ©-([I]) yield the desired alge- 
braic theory for the leading-order correction to go (r) in 
the particular single square-well case. Notwithstanding 
its simplicity, the resulting theory is in striking quantita- 
tive agreement with "exact" Monte Carlo results (Fig.[T]). 

The above model can be directly extended to encom- 
pass multiple steps of arbitrary widths; thus, for a po- 
tential with M steps, 



gi (r) = -figo(r) 



(5) 



where — £j is the energy of the ith step, and v^'(r) is the 
volume of the intruder particle inside the ith step region 
(defined analogously to v*(r) in Fig.[3j. A representative 
application of this model for M = 2 is shown in Fig. [2] 
With the exception of the more pronounced deviation 
near the contact distance r/a w 1 (chiefly due to the 
solvation approximation, and not to the model itself), 
the volume-exclusion theory is again in good quantitative 
agreement with direct simulation results. 

Since corrections at a given order to the structure fur- 
nish corrections at next order to thermodynamics, it is 
of interest to investigate the consequences of the above 
first-order structural theory to second-order corrections 
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FIG. 4: Relative error (A" 2 oLv - A e 2 xact ) /\A e 2 xact \ in the 
second-order correction to the free energy, viz. A 2 /N = 
(p/4) J dr it(r) <;i(r) [3], versus density, where the super- 
scripts specify the use of Eq. ([Taj) [A e 2 xact ] or Eq. IgaJ [A s 2 olv ] 
for the computation of gi (r). From top to bottom: Discrete 
core-softened (CS) model of Figure [2] and single-well poten- 
tials with ranges r* /a — 1.7, 1.5, and 1.3. 



in the free energy, Ai [14} . For the models studied, 



the theory predicts A^ with relative errors of the order 
of 10% at densities per 3 > 0.8, the source of discrep- 
ancy at lower densities being traced back to the under- 
lying solvation approximation (Fig. 3]), as already an- 
ticipated. It should be emphasized, however, that the 
volume-exclusion model above is in principle justified 
only at high densities, and hence the gradual breakdown 
of the solvation approximation at lower densities is con- 
sistent with the scope of the algebraic theory. 

To this author's knowledge, the foregoing is the first 
structural perturbation theory of liquids to offer analytic 
corrections to go(r) of such conceptual and algebraic sim- 
plicity that nonetheless yields meaningful results at high 
densities (compare with the superposition approximation 
14j; for algebraic corrections to thermodynamic proper- 
ties at first-order, where only go(r) is required, see the re- 
cent work of Ben-Amotz and Stell (23[). This is achieved 
through the introduction of two novel ideas: the solva- 
tion approximation (Eq. @) - in principle applicable to 
any dense liquid with either discrete or continuous poten- 
tials - and a simple volume-exclusion model for three- 
particle correlations in liquids with discrete potentials 
(Fig. [3]). Given its ability to easily handle wells with 
multiple depths, it is hoped that the methodological de- 
velopment put forward in this paper will help elucidate 
the role of length scales in liquid models with thermody- 
namic anomalies 0, [Io| , as well as in perturbative treat- 
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